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ABSTRACT 

A method, using simxxLtaneous equation solving technique, 
has been developed to solve multicomponent distillation problems 
and then extended to simulate a crude-oil distillation column. 
The crude oil is assumed to be composed of a nmber of pseudo- 
components whose thermodynamic properties are correlated with 
their average boiling points and specific gravities. In this 
method, the equilibrium, overall material balanpe and enthalpy 
balance equations have been solved simultancoiasly with L, V, T 
as independent variables using Broyden's procedure. Component 
material balance equations were solved using modified Thomas 
algorithm as have been proposed by Hoefling and Reader, 
Fractionation problems having large number of components of 
wide volatility range can be solved by this method. Its 
convergence characteristics have been found to be satisfactory. 
The memory requirement is also low which is important for 
columns with large number of trays and feeds containing many 
components. 



CHAPTER 1 


INTRODUCTION 

Petroleum industry is one o£ the biggest among the 
chemical industries with an annual throughput of about 40 
million tonnes of crud© oil •The entire crude oil to be 
processed by a petroleum refinery has to pass through the 
crude distillation unit or more commonly called toping 
unit which makes it the most important of all process equipment 
in refinery. Therefore the importance of rating or design 
of distillation columns processing crude oil or any other 
petroleum fractions is imaginable* Before the advent of 
modern fast computing machines, people used empirical methods and 
rules of thumb to design crude distillation columns. With tiie 
rapid growth of electronic computers, a numberof methods of 
multicomponent distillation calculations have been published 
(Amundson and Pontinen, 1958, gujata, 1961, Rose et al* 1958, 
Tomich, 1970, Naphthali and Sandholm, 1971 etc.). Consequontly 
interest to simulate or design distillation columns processing 
petroleum fractions using computers has also grown. But till 
date very few papers have been published in this regard 
compared to the steady supply of p^ers on multicomponent 
distillation in Chemical Engineering journals* The refineries 
still continue to use empirical correlations and thumb ruiles 
for design of their distillation units. It is not the lack 
of interest in use of modern methods of computation on the 
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part of petroleum refiners which has prevented them from using 
these tools, rather^, it is the lack of appropriate procedures 
which are ^plicable to their problems. Although crude 
distillation like most other distillation problems is a 
multicomponent, multistage process, yet it differs significantly 
from other conventional columns* Whereas, in conventional 
multicomponent distillation units, the feed consists of a 
specified number of identifiable chemical compounds to which 
thermodynamic properties can be assigned uniquely, it is neithei" 
practical nor possible to identify the constituents of a crude 
oil. Therefore special procedures are required to characterize 
the crude oil to enable the estimation of thermodynamic 
properties s\x:h as vapor-liquid equilibria, enthalpies etc. 

In addition to the above difficulties the complex nature of the 
tower configuration having side strippers, unconventional type 
of reflux, heat removal at ihtermediate plate create further 
conputational problems, wide boiling range of the petroleum 
fractions to be distilled from crude oil may lead to additional 
ptoblems in convergence. 

In the present study, attempt has been made to model a 
petroleum fractionator by considering the feed to be composed 
of a finite number of pseudo-components to which thermodynamic 
properties can be assigned on the basis of their average 
boiling point and density. Once the feed is characterized in 
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this fashion^ any multicomponent calculation procedure 
could be used for rating of a refinery distillation colunn. 

In this study> a modified version of Tomich (1970) method 
has been used# 

Ihe whole thesis has been organised in the following 
way, ihe second chapter contains* literature reviev/ of 
computational methods for multicomponent fractionation in 
general, Ihe third chapter describes features and problems 
specific to crude oil distillation. The formulation of the 
model by the present method is included in chapter four. 'Ihe 
fifth chapter contains the results of computation obtained by 
the present model and compares them with those obtained with 
a standard commercial package,' The last chapter presents the 
conclusion drawn and recommendations for future work. 



CHAPTER 2 


LITERATURE REVIEW 

In this chapter major developments in the field of 
simulation of multicomponent fractionation have been reviewed 
with an emphasis on Newton-Raphson technique used in solving 
these problems. 

2,1 Simulation of Multicomponent Fractionation : 

Since the advent of modern fast computers there have 
been many articles published in literature which use plate to 
plate rigorous calculation. This trend in multicomponent 
calculation started with the historic article by Thiele and 
Geddes (1933), This study used the sequential method of 
calculation where starting from both ends of the coliamn plate 
by plate calculations were made sequentially until the feed 
plate was reached. The convergence in these type of plate to 
plate calculation methods was usually slow which was promoted by 
a modification introduced by Lyster et al. (1959) called the 
*e method* of convergence promotion. Pondon et al* (1973) 
proposed a multi—© method where a different theta was used for 
each plate. The application of © - method of convergence has 
been demonstrated by Holland (1981). Things turned in a 
different direction using more elegant mathomatical techniques, 
when Amundson and Pontinen (1958) proposed a formulation in 
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which one writes the necessary material balance equations in 
a compact matrix form which can be solved by efficient 
algorithms such as the Thomas algorithm* A variation of this 
type of method was proposed by Sujata (1961) which is popularly 
known as the *sum rate* or SR method/ whereas the original 
one is termed as the ‘Bubble point* or BP method* These two 
methods differ in their way of updating temperature in each 
iteration* The BP method utilises bubble point calculation 
and computes vapor flow rate from enthalpy balance equations/ 
whereas the ’SR* method ccsnputes the total vapor flow rate by 
summing up the component flow rates vhich are obtained by solving 
tridiagonal coefficient matrix and then evaluates the temperature 
profile from enthalpy balance equations. The *BP* method is 
generally applicable to close boiling mixtures/ hence best 
suited for distillation problems, whereas the 'SR* method can 
solve wide boiling mixtures/ like the kind faced in absorption 
problems* The applicability of these methods has been discussed 
by King (1971). Rose et al* (1958) introduced a different 
approach in multicompxDnent calculations popularly known as 
'Relaxation technique*. It follows the transient behaviouir of 
separation processes approaching steady state and uses a 
relaxation factor* Ihe successful application of this method 
depends upxsn the choice of optimal value of this factor* 

Tomich (1970) and N aphtha! i-Sandholm (1971) proposed tvjo 
methods which used multivariate Newton-Raphson (NR) technique 
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in solving multicomponent problems. These two methods can be 
termed as simultaneous methods where material and energy 
balance equations along with the equilibrium relations are 
solved simultaneously. This is in contrast to the ‘BP* and ’SR’ 
methods which are successive iteration methods. The simultaneous 
methods are considered to be superior to successive ones since 
in the former, all the necessary variables in the functions 
are corrected simultaneously as distinguished from the latter 
class of methods where only one set is corrected at a time. 

These two methods need be described in greater detail as our 
present method owes its origin to these two methods. These have 
been presented in Sections 2,2 and 2,3, Ketchum (1975) proposed 
a new method which combines the salient features of both NR 
technique and Relaxation method. This method strives to attain 
the speed of convergence of NR method while retaining the 
stability of relaxation method, HlavaceJc (1981) introduced one 
parameter imbedding technique in multicomponent calculation. 

Also different authors have tried to use NR method by grouping 
the variables and error functions in different ways. Notable 
among these is a paper by Perraris et al, (1982) which is 
considered to be suitable for systems having one component of 
either very high or very low volatility. This method can be 
applicable to systems having steam present as one component 
such as steam strippers. This method can be considered as a 
variationof Naphthali-Sandholm technique, Chakraborty (1982) 
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has shown several different variants of Naphthali-Sandholm 
method. 

2.2 Tcynich*s Method ; 

This procedure basically solves component material 
balance equations with a tridiagonal matrix algorithm. Heat 
balance and summation equations are handled with a quasi - 
Newton method, A notable feature of this method is that^ in 
a mathematical sense, all equations are solved simultaneously, 

2.2.1 Description of the Model Column: 

The model column for mathematical simulation is shown 
in Figure 1, Because of the general form, the model represents 
all types of stage-by-stage separation processes, though in 
the figure a distillation column is shown. Feed streams can 
be introduced on each tray of the model. Vapor and liquid 
stream can also be withdrawn from each plate, similarly heat 
can be added to or withdrawn from each tray to simulate 
reboilers, side heaters, interchillers and condensers, 

2.2.2 Model Equations: 

The basic equations for column simulation are derived 
by making material and heat balances around jth tray of the 
model (Figure 2), Heat and mass balances plus equilibrium 
relationships comprise five sets of equations to be solved. 

(i) Overall material balance equations: 

L. , - (Its .) L. - (1-tfi ) V, + V. + F. = 0 
3 3 3 3 3+1 J 



feed 



M 

f 

N-1 


FI6.1 
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- L. 
J+1 j 


V . = 0 
J 


j = 1 

j = N 


C2-.1) 


L. T - V. - L. = O ; 

J J 

where L = liquid flow rate 
V = vapor flow rate 

W. 

Sj *= fraction of liquid drawn off = ^ 

3 L. 

U. 

S. = fraction of vagpor drawn off = ^ 

^ J 

F = feed rate 

N * total number of plates 

Si±>script j-ly 3 , j+l represent the (j-l)th, jth and 
(j+l)th plate respectively, 

(ii) Vapor-liquid equilibrium relations: 

° ; 1 4 i ^ M 

and 1 ^ j 4 ^ N (2.2) 

where K = distribution coefficient 

y = mole fraction in vapor phase 

X = mole fraction in liquid phase 

M = number of ccanponents 

Subscript i refers to the ith component* 

(iii) Component material balance equations: 


y< .1 4 + 1^4 1 ^4 1 ■! + ^4 Z. . 

j "fl j+l/i J— 1 J— 1/ i J J/i 




- (l+Sj)Vj = o 


1 i ^ M 
and 2 ^ j ^N-1 
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where 


^3+1,1 - ^3 ^3,1 - ''j “ ° ' 

and 


L . T X . . . - L . X . . - V. y . , = O ; 
0"“1 J""!/ i J 3/i J 3^^ 

and 


1 ^ i 
j = 1 

1 4 1 ^ M 


(2.3) 


Z = mole fraction in feed 

(iv) Summation equations: 

M M 


‘j' jli |i 


(2.4) 


(v) Heat balance equations : 


= rv. ,1 + L . - h . . + F . h 

I j +1 j +1 3-1 3-1 3 


-(14«,)L. h. - (14S )V. h1- 

3 3 3 3 3 3 


Qj - 0 


I 2 J N-1 


-3 = ''3+1 ”j+l - - '’fi - “c “ ° 

; 3 »= 1 


E. = L. , h. . - L.h. - V.H. + Q = O 
j 3-1 3-1 3 j 3 3^ 

3 = N (2 .5) 

h * liquid stream enthalpy 

H =* vapor stream enthalpy 
F 

h “feed stream enthalpy 
Q » heat removal rate at intermediate plate 
* condenser heat duty 
=® reboiler heat duty 


where 
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2 ,2 *3 Method of Solution: 

Equation sot (2,2) can be substituted into set (2,3) 

to eliminate x. j and reduce the component material balance 
3 / * 

equations to the tridiagonal matrix form for the case of simple 
distillation columns* 


'b^ 1 


'^l,i 


'°l' 

^2 ®2 ^ 


'^2,1 


°2 

C, 1 



ss 

D-> 

3 3 

mm mm mm 


3, X 


3 

S-l 1 


V 

N-1, X 

1 

1 

j 

%-l 



\,± 




. 1 i ^ M 

( 2 . 6 ) 

In arriving at the above equations, use has been made 

of the fact that v. ^ = V. y . . and 1 . . = L. x. . where v. . 

3#i 3 3/1- 3 3#'i 3/i 

and 1 . j are the vapor flow rate and liquid flow rate of ith 
3/ * 

component leaving jth stage respectively. 

Also = - (R+1) for total condenser, 

where R is reflux ratio • 

Aj,,= LjAVj 
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lliis set of equations (2.6) is a linear or nearly linear 


set when K-values are composition independent or weak-functions 
of composition. The procedure to solve all the equations is 
summarized stepwise as follows: 

1. Initial values of Tj and are assumed for 
all 3 • 


2. Lj s are computed fxx>m equation set (2,1). 

3. Equation set (2,6) is solved using Thomas algorithm 
(Lapidus, 1962) or any other efficient algorithm to 

give V. j's and hence y. j ‘s. 

J/i J/3. 

4, Equation set (2,2) is solved to give x. j ‘ s* 

3/ i 

5, These results are substituted in equation sets (2,4) 
and (2,5) \4iich are solved simultaneously by Broyden*s 
procedure to give a new set of T^ and V^, 

6, With these new values of T. and V., calculations 

J 3 

from step 2 through step 5 are repeated until 
convergence is achieved. 

Convergence criterion used in this case is 


C2.7) 

j=l 

where £ is preassigned tolerance limit, 

2,2,4 Critique of Tomich*s Method: 

It is a general method which can be used for any type 
of multicomponent fractionation problem such as distillation^ 
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absorption etc* By using simultaneous equation solving 
procedure, it eliminates numerical unstability which occurs 
in sequential approach due to build-up of truncation errors. 

Use of quasi-Newton method reduces computation time. Since, 
for each iteration, Jacobian and its inverse need not be 
calculated. lOmich^s method requires less memory space 
compared to other methods using NR technique like Naphthali- 
Sandholm*s method because of less number of iterated variables. 
However, this method fail to converge for highly nonideal 
sol utions , 

2»3 N aph thal i— Sandholm Method ; 

In this procedure, the equations of conservation of 

mass and energy and equilibrium relations are grouped together 

stagewise and then solved by linearizing them. Ihe resulting 

set of equations has a block tri-diagonal structure. 

Figiare 2 represents a model plate whose side streams 

are specified by the ratio of the side stream to the stream which 

remains after they are withdrawn. In this formulation 1. j# 

V , T, are the independent variables vhereas L V. representing 
3 3 3 

/ 

the total phase flows within the column are dependant variables. 
The plates may have a feed stream introduced to it and there 
may be heat removal from or heat addition to it as before. 
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2,3.1 Model Equations: 

(i) Enthalpy balance equations; 




(1 + S.) V. , H. . + (1 + s.) 1. . h. . 
3 3/3. 3 3/3^ j/3. 


^j-l/i ^j-l#i ^j+1, i ^j+l,i “ ^j,i ^j,i3 


1 ^ i 4 M 
and 2 ^ j ;< N-1 


( 2 . 8 ) 


where H *= enthalpy of components in vapor phase 

h » enthalpy of components in liquid phase 
P 

h « enthalpy of components in feed stream 
f = component flow rate in feed stream 

Sj - fraction of vapor withdrawn from jth plate 

= u / V 
3 3 

Sj = fraction of liquid wL thdrawn from jth plate 
= 

(ii) Component material balance equations; 

M. . = (1 + S ) V . + (1 + s ) 1. . - V 
3/3. j 3,1 3 j,i 3+1/1 


j-1,1 j,i 


1 ^ i ^ M 
and 2 j ^ N-1 
(2.9) 


(iii) Vapor-liquid equilibrium relations; 


Ej j K. , V. 1 . . 


(1-E, ,) v,_,, ^ V, 


- V. + 


3+1 


/* 1 ^ i ^ M & 2 ^$C 3’ N-1 


( 2 . 10 ) 
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where E *= Mxjrphree plate efficiency 

Q. . is derived from the definitions of the vapor phase 
J/ ■*- 

murphree plate efficiency, that is 


E . = 

k, , X,, » y 


J/3. 3,i 


j+l/i 


( 2 . 11 ) 


or 


E . . K . . X . . 

3,i 




0 


( 2 . 12 ) 


where x and y represent mole fractions in liquid and vapor 
phase respectively. In terms of component flow rates it becomes 


V 


'j/i 


K 


, . --Iti - -jiJ: + ( 1 - E . . ) — iiiii 


1/ i Lj 


1 / 1 ' 


V 


O 


j+1 


(2.13) 


Mioltiplying both sides by equation (2.10) is arrived- 
2.3,2 Method of Solution: 

Ihese functions apply to all intermediate plates of the 
column. Ihey can also be applied to partial condensers and 
reboilers provided reboiler heat duty and condenser heat duty 
are known. There can also be different other specifications 
possible \hich are detailed in the paper (Naphthali and 
Sandholm, 1971). There are (2 m+ 1) number of variables and 
same number of functions for each stage or a total of n(2m+ 1) 
variatie and equations. These are grouped stagewise. In 


compact form theS'e can be written as 
F (X> = O 


( 2 . 14 ) 



17 


where 


^ “C^l' ^2' ^3' ' 

F =[?!, Fj, F 3 F^J ^ 

•m T 

^■i ^r^-i 1' 1. ,, 1. 1 

J L 3^2 2 ' j,l' j ,2 j,Mj 

F. = Te / M , M. Q. w Q. ,..^Q. "j 

J L J J/1 J/2 3,M' 3/1' 3/2' 3/MJ 


-Tnus using this notation the NR method becomes/ 

2"''^ = 5 "> + A X 


wh ere 




= 5 "> + A X 

= -5 -1 i ■" 


(2.15) 

(2.16) 


^ ^ ibh- 1 . . , . — m 

X xs the correction added to X to get the new values 

•M M j 

of the variables x • Here superscript m and m+1 refer to 
iteration nimber. If the functions are linear/ this correction 
would make the value of each function equal to zero. Since 

•w m I 1 

the functions are nonlinear/ A. X is just an approximation 
to the needed correction. J is the Jacobian matrix that is the 


matrix of partial derivatives of all the functions with 


respect to all of the variables at present value of the 
variable X Hence 


^3_ 

3 Xi 3 X 2 

= _ _ ^^2 ^^2 ^^2 

BX2 

. ♦ 

m 

^^1 


(2,17) 
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where 



1 




i, 2 


^ x .^1 




iy 1 


B F. T 

^ ^j,2M+l 


2m +1 ■ 2M+1 

0 ^j,l ^^j,2M+l 

( 2 . 18 ) 


It is netted that the functions for plate j involve only 
variables of j“l, j and j+1. Thus the partial derivatives of the 
functions on jth plate with respect to the variables on all 
plates other than these three are zero. So the Jacobian matrix 
takes the structure of a block tridiagonal form 


Q 


J 





c . 

J 



^-1 ®N-1 S-l| 

L ^ J 

(2.19) 

where A, B, C are individually matrices of the order (2 m+1). 


a. = ( 1 — 


); 


B 


©F 


C , = ( 1— 


) 
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Detailed structure of these siib-matrices can be found 
in the original paper referred above* 

Newton Raphson method by itself is often unstable, 
specially when the initial guesses are far from the solution* 

To reduce the chances of divergence, a damping method may be 
used in which, 

- m+1 ^ - m ^ ^^-m+1 (2*20) 

The value of the damping parameter, , should be chosen 

—mtl 

in such a way that X. minimized. For this, any one 

dimensional -search technique like Golden Section method can be 
used* Even with ^ -modification it is important that initial 
guess is not far from the solution. One procedixre to generate 
a reasonable initial guess is to assume, a constant molal overflow 
and a linear temperature profile* Then, with composition 
independent K-values, solving the component material balance 
equations, component liquid and vapor flow rates through the 
column can be generated* 

2*3*3 Crique of Naphthali-Sandholm Method; 

This is a very general method which can be applied to 
any type of multicomponent calculations without having to 
modify it. Volatility range does not affect the method either. 
Without requiring any modifications, highly non- ideal solutia:is 
can also be treated because the presence of non- ideality is 
accounted for rigoirously* various types of specifications 
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can be taken into account in a rigorous manner through the 
equilibrium relations. Feeds can be added to or side streams 
can be withdrawn from any plate* Also heat can be added to or 
can be removed from any plate. The method is based on a 
linearization of distillation equations. So convergence 
accelerates as the solution is approached. No computational 
difficulties or complications arise from component flow being 
small. The primary disadvantage of this method is the huge 
amount of storage required. The size of the Jacobian matrix 
in this method is N(2 m+1) xN(2m+ 1) . Though with efficient 
m«nory saving algorithm like block Thomas algorithm 0 y \\^ 

Cj sub matrices need to be stored^ with interconnected 
colmns, the storage requirement goes up. 

2.4 Interconnected Column : 

The simulation of interlinked colunns can be done by 
the same algorithm treating all the columns present as a single 

column. Hbwever, the statement “Functions of a plate contains 

A 

only the variables of that plate and the two adjacent plates'* - 
does not hold good- any more due to the interconnections. 

Because of this, the Jacobian matrix contains some off-diagonal 
elements besides the tridiagonal band affecting the sparsity 
of the matrix* Stadtherr (1979) proposed a technique to 
maintain sparsity in process design calculations. Hidalgo et al* 
(1980) proposed an optimal arrangement of simultaneous linearised 
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equations for general systems of interlinked multistaged 
separators. They developed an algorithm which automatically 
arranges the linearized equations so that for a Newton -Raphson 
type of solution procedure (for example^ Naphthali-Sandholm 
technique), a minimum or nearly minimian number of non-zero 
blocks occur outside the tridiagonal band. In this way, 
computationsd memory and time requirements are minimized. 
Hoefling and Seader (1978) proposed an algorithm to 
solve the tridiagonal matrix having offdiagonal elements, which 
is a modification' of original Thomas algorithm, Seader et al. 
(1980) proposed an optimal arrangement of equations for inter- 
linked columns to minimize the non-zero elements, Stadtherr 
et al, (1982) looked into various solution techniques for 
sparse matrices, common in multistage separators, Stadtherr 
(1982) in another paper proposed an algorithm which rearranged 
the block tridiagonal matrix having off-diagonal elements into 
bordered block tridiagonal form and then applied block 
triangular splitting technique to solve it. This algorithm 
is comparable to Hoef ling-Seader *s algorithm in efficiency. 



CHAPTER 3 


CRUDE OIL DISTILLATION 


Crude oil distillation is one type of distillation 
which still needs special attention/ inspite of many general 
purpose package available commercially because of complexity 
of crude composition/ unconventional type of column con£igU“ 
ration and its operation. In this chapter the process of 
crude oil distillation will be described briefly and then 
the major methods available for its calculation will be 
discussed, 

3.1 Distillation Column : 

Unlike the mixtures treated in usual distillation 
columns which are made up of a few discrete components/ crude 
oil consists of a very large number of hydrocarbons all the 
way from methane to materials having 70 or more carbon atomS/ 
the exact composition of which is nearly impossible to determine 
Therefore the way to characterize crude oil or a petroleum 
fraction for the purpose of its fractionation is to determine 
the boiling point curves from results obtained in laboratory 
test distillation units. Depending upon the operating condition 
there can be different types of boiling curves like TBP/ ASTM/ 
EFV. The specifications and operating conditions of the 
apparatus to obtain these curves are given in detail by Nelson 
(1958), True Boiling Point or TBP distillation is usually 
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Condenser 



F16.3 


A TYPICAL ATMOSPHERIC CRUDE DISTILLATION TOWER 
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run on crude oil and is indicative of its composition, whereas 
ASTM is run on different petrol emi fractions and used by 
designers to specify a certain fraction# EFV distillation is 
seldom run due to its high cost* All other properties are 
assigned to the fractions by correlating these properties 
with the average boiling point and density of the fractions. 

In crude»oil distillation, the crude is preheated in 
a furnace and then flashed in the flash-zone of the distillation 
tower. A schematic diagram of a typical refinery crude 
distillation unit is shown in Figure 3. A small amount of extra 
vaporization^ called overflash, is employed to provide adequate 
reflux between the flash Zone i.e* the point at which partially 
vaporized feed is introduced and lowest side stream product 
draw tray. All the heat supplied to the column comes through 
the preheated feed, as there is no reboiler in the tower. Open 
steam is introduced into the column to strip the lighter 
components from the bottom prodxxrt as well as to bring down the 
operating pressure effectively. Liquid products are withdrawn 
from different points of the column through side streams. 

Lighter components are stripped by open steam from these side- 
streams in side-»strippers and then sent back to the column# 

Apart from top tray reflux, heat removal at intermediate points 
is attained by withdrawing an internal liquid stream from the 
tower, cooling it and then returning it to the column# The 
lightest distillate is not always condensible at the conditions 
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of temperature and pressure in the reflux drum and thus may 
be treated as two distillate fractions, one vapor and the other 
liquid# 

3.2 Methods of Calculation ; 

'There are two major trends in crude-oil distillation 
calculation* The first one is empirical in nature and is 
based upon Packin' s method (1941) and second one utilizes tiie 
pseudo-component concept. We briefly discuss these methods 
below. 


3*2,1 Bnpirical Methods: 

In this method, for petroleum fractions, two terms are 
used to describe product composition and the degree of separation 
between two adjacent streams, ASTM boiling range defines the 
general composition of the fraction and is usually one of the key 
specifications for most distillates. The second term (5-95) gap 
signifies the relative degree of separation between adjacent 
streams. It is defined as the difference between 5 volume 
percent ASTM temperature of the heavy fraction and 95 volume 
percent ASTM temperature of the adjacent light fractions. 


(5-95) gap = (t^ - t^^ ) ASTM, degrees FarenhQit 

H L 


(3#1) 


In this method, the degree of difficulty of separation is 
defined as the difference between the 50 volxome percent 
tenperatures of the fractions under consideration. The ‘separation 
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capab^ility • of the systan is defined as the product of the 
reflux to feed ratio at the uppcar draw tray (as calculated on a 
volumetric basis) and the number of trays in that section. 

The product is designated as F-factor. In the systems where 
pump— around reflux is used for heat ranoval the trays employed 
in this service are regarded as one- third of an actual 
fractionating tray. The degree of separation attainable and 
the (5-95) gap are defined as functions of the separation 
capability of the system (i,e. P-f actor) with parameters of degree 
of difficulty of separation (i,e, ^^50 ASIM) . The functions 
Or correlations are basically empirical in nature. Later Packie*s 
method has been extended by Watkins (1969), 

3,2,2 Pseudo— Oamponent Methods: 

In this method/ the crude mixture is considered to be 
composed of an umber of pseudo-components or fractions/ each 
of which having been assigned average properties such as boiling 
point/ molecular weight, density, critical and other properties. 
The basic steps in this method are briefly stated below. 

From the TBP-curve of the crude oil feed (see Figure 4) 
a number of constant temperature cuts are stepped off so that 
areas above and below the line between it a nd TBP-c;irve are 
equal (for example area abc = area ede in Figure 4), These 
cuts are termed as pseudo-components. The present components 
are considered to have average boiling points equal to those 
constant temperatures. For step abede (Figure 4) the pseudo- 
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component has the boiling point T^. Tlie number of these cuts 
is decided upon the preciseness of the separation desired. 

Hie width of these steps need not be equal. In the regions 
where the TBP-curve is sharp^ the step width can be small 
whereas in flatter regions of the curve steps can be wider 
i.e, number of pseudo-components can be fewer. The density- 
midpercent curve is drawn from the experimental data on the 
cuts and hence each cut is assigned an average value of density* 
The molecular weight, watson characterization factor, critical 
parameters, acentricity factors of the cuts or the pseudo- 
components are determined from various correlation involving 
the boiling point and density. 

To get the mole fractions of each-component in the 
feed mixture, the volume of each cut is converted to weight. 

The weight of each cut is then divided by its molecular weight 
to get the moles of each cut and they are summed up to give 
total molar flov; of the feed. Then individual moles of each 
component is divided by total number of moles in the feed to 
give its mole fraction in the feed. Thus having got the mole 
fractions of the pseudo-components in the feed and having 
assigned molecular properties to each pseudo-ccanponent, any 
short cut method of multicomponent calculation or any rigorous 
technique as described in Chapter 2 can be applied for 
simulation of crude oil distillation. Prater and Boyd (1955) 
used this pseudo-component concept along with short cut methods 
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to design crude oil distillation twoer, Cecchetti et al,(1963) 
first solved a crude oil distillation problem using this 
pseudo-component concept and a rigorous technique of multi- 
component distillation. They used Thiele-Geddes calculational 
procedure with 9 - method of convergence. They split the 
crude oil into 34 pseudo-components. In their treatment an 
analogue of the actual column was fourd by trial. Analogue 
of the actual column is defined as a column having s-uch number 
of theoretical plates that gives the same product stream as 
obtained from the actual column. Later, Hess et al, (1978) solved 
the same problem using "2 n Newton Raphson technique". The 
difference in their approach was that whereas Cecchetti et al, 
assumed water vapor to be present only in the vapor phase 
throughout the column, Hess et al, considered water to be 
present in both phases, 

3*2.3 Other Methods: 

Besides above two approaches, a few others have tried 
somewhat different methods to solve crude distillation problems. 
Notable among these is the work of Nikola j (1970) who solved 
crude oil distillation problem by making different modioles for 
different sections of the column, Ihese modules can be solved 
independently* Material and energy balance calculations were 
done around each module. The properties of the product streams 
of the modules were computed by using empirical correlations. 
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By making different combinations of these modules, different 
problems having different tower configurations can be solved. 

A Critique on Crude Oil Distillation Calculations : 

It should be anphasized that no method, including the 
pseudo-component method, is exact because of the assumptions 
necessary to assign average properties to each of the pseudo- 
components* In principle, the pseudo component method will 
approach exactness as the number of cuts approaches infinity. 

It is, however, not practical to make more than about 50 pseucSo- 
component which may not provide a reasonable approximation to 
the actual feed. To assign thermodynamic properties to these 
components, one uses generalized correlations which may 
introduce further uncertainties in calculations, since it is 
very difficult to acc irately correlate these properties with 
the available specifications. Precise vapor-liquid equilibrium 
data are very difficult if not impossible to obtain. It is 
nearly impossible to define and thus evaluate stage efficiencies. 
This problem can be tackled in a round about way by defining 
an analogue column for a particular distillation unit. However, 
once an analogue of a particular column is found by trial 
and error it can be used to solve any problem with that column. 



CHAPTER 4 


MODELLING OF CRUDE OIL ETSTILLATION 

Naphthali-Sandholm is by far the best method available 
for multiconponent calculations. But its use is limited 
because of its large memory requirement. The Jacobian matrix 
in this type of formulation is of the order of n(2m+ 1) and 
for systems having a large number of components it can create 
a great deal of problem. For simple distillation problems 
where there are no interconnections present, the Jacobian 
matrix will have a regular shape of block tri-diagonal form 
because of stage-wise grouping of variables, Edfferent 
efficient memory saving algorithms can take care of this type 
of problem without much difficulty. But with interconnections 
present, as in case of a colmn with side strippers, elements 
also appear in positions other than the tri-diagonal band in 
the Jacobian matrix resulting in its greater fill-in. In crude 
oil distillation due to intermediate refluxes, the structure 
of the blocks also get changed. For these reasons, Naphthali- 
Sandholm method can not be applied very easily to solve this 
type of problem. 

Another alternative is Tomich*s method which has low 
memory requirements. Its Jacobian matrix is of the order of 
2N but without the tridiagonal band structure. As mentioned 
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QaiTlier (Section 2*2,4)/ its convergence characteristics can be 
poor, Chajcraborty (1983) has reported the use of some variants 
of Naphthali~Sandholm method which have better convergence 
characteristics than TOmich's method^ thou^ poorer in convergence 
with respect to* the original Naphthali-Sandholm method. In 
the present work, we have chosen the functions similar to 
those suggested by Chakraborty. The dependent variables/ 
however/ no longer depend only on adjacent stages which leads to 
the loss of block tridiagonal matrix characteristic of Naphthali- 
Sandholm method and its variants. The present formulation 
deviates from that of Tomich in the sense that while Tomich 
assumed and corrected temperature and total vapor flow rate on 
each stage simultaneously and subsequently calculated liquid 

flow rates, in the present case, all the three variables T.# V. 

J j 

and Lj are assumed and corrected simultaneously, 

Althou^ in the present study, the formulation as 
discussed above has been used, advantages of Naphthali-Sandholm 
method cannot be discounted. Ihe non-availability of enough 
storage in computer prevented the testing of Naphthali-Sandholm 
method for crude-oil distillation problem in the present study. 
But the formulation has been included in Appendix A and it is 
hoped that this approach will be tested at a later date. 
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Modelling o£ a Distillation Column VJithout Interconnections ; 

Before attempting to model crude oil distillation 
unit/ the present approach has been applied to a conventional 
multicomponent distillation column without interconnections. 
Since such colmns can as easily be simulated by Naphthali™ 
Sandholm or other methods, this exercise enables us to establish 
the validity of this approach to at least such problems. 


4*1,1 Ihe Model s 

The model plate of figure (2). is considered once again. 

The equations describing the plate were presented as equatios 2,1 

through 2,5, Substituting the value of 1 . ./l. from equation 

1/ ^ J 

(2,2) into equation (2.4), the following equation is obtained: 


1 

V. 

3 


M 


21 

i«=l 


* 0 


j l/2/«*. 

(4.1) 


Following three functions are chosen for iterative 
calculations. 

The equilibrium relations: 


^3 ‘ 


^3 


M 


i=l 


'K 


- 1 ) v. 


j/ i 




j “ 1 / 2 / • . » / N 


(4.2) 


For partial condenser and partial reboiler the equilibrium 
relations remain unchanged. 

The overall material balance equations: 

- (1 + = 3 ) - (1 + Sj> Vj + ; J = 2,3,. -H-l 
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M 


j 


M 


j 





/• 3 

/• 3 


« 1 
W N 


(4.3) 


The enthalpy balance ecjuaticwis: 

M 

V = (1. - 1 4 ~ (1-te.) 1. , h. . - (l+S.) V. . H. . 

J ^ J j/i' j/i j J/ i 


1*^1 TT 

+ V. ^ . H 4- f . h . . ) 

3 “^ 1 / 3 . 3 /^ 3 /^ 


j ~ '2. f 3/ • . ./N“l 


The enthalpy balance equations for j=l are replaced by 
the following specification equations, 

- RV^ 

or Ej^ = V. - D 

where R = reflujc ratio 

D = distillate flow rate 

For the reboiler (j=^^) specification equations become, 

% ' '^N ' 

or (4.4) 

where T* and L* are the specified temperature of and liquid 
draw off from the reboiler respectively. There can be other 
specifications also possible, like, reboiler heat duty or 
reboil ratio. For total condenser, the equilibrium equation 
is replaced by a specification equation. The specification 
in this case may be condenser temperature or condenser heat 
duty, 

4,1,2 Method of Solution: 

For the three functions (Equations 4,2, 4,3 and 4^,4) 
at each plate, the chosen independent or iterated variables are 
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T., L. and V., 
J j J 


where 


In function form these equations become, 
F (x) = 0 


and 


F » 1 


^2' 

F .** 

3 

[Sj, 


.. 



X — 

[=^1' 

^^2' • * * * • 

X. = 



3 


- 




For a given set of values of the independent variables, 

the dependent variables v. . and 1. . can be found in the 

J/ i J/ i 

following way* The vapor phase component flow rate v. . can 

1/ i 

be computed by solving equation (2,6) p 1. can then be 

3f ^ 

calculated using following relations 


1 . , = A. . V . . 

j, i 3,3- 

where A . . = L ./v . K . . 

3# 3- 3 3 3/3. 


(4.5) 


There are three equations at each plate, namely, equations 4.2, 
4.3 and 4.4 and also there are three unknowns, namely, Tj, 
and Vj for which those equations are to be solved. For a total 
of N plates, therefore, there is a set of 3M equations in as 
many variables to be solved simultaneously. The major steps 
involved in the solution are as follows: 

Step 1 ; Assume initial values of T^, Lj and V^s which form 
the group vector X 

Step 2t With these values of T^., L^. and V^s solve equations 
(2,6) and (4.5) to generate component flow rates 
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throughout the column. Thomas algorithm is used 
to solve equation set (2,6), since, the coefficient 
matrix has the tridiagonal band structure. 
step 3; Equations (4,2), (4.3) and (4,4) can now be solved 
using Broyden's algoritlwn which is discussed next, 
Broyden’s Method: 

Broyden’s method is based on the numerical approximation 
of the partial derivatives appearing in the Jacobian matrix 
(Broyden, 1965), This approach permits the inverse of the 
Jacobian matrix to be updated after each trial. Thus it is 
necessary to compute partial derivatives in the Jacobian 
matrix only once. As approximate values for the partial 
derivatives are used, it requires more iterations to converge 
to the solution. However, since it avoids computation of 
partial derivatives and its inversion with every iteration, it 
saves a lot of computer time per iteration, specially when 
numerical differentiation is used for the calculation of 
Jacobian matrix. In this method, the following calculations 
are performed to solve the set F (x) =0, 

S., M. and E- which are the elements of F. are really 
3 3 3 3 

the errors or residuals in the equilibrium relation, material 
and heat balance equations and are equal to zero only v/hen the 
solution is reached. These functions are expan ded in Taylor 
series and terms containing derivatives of higher than first 
order are dropped to get 
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J Z5v X = - p 


(4.6) 


where j = 


5^1 


'c) ^ 


ifi' 

^ ^ i 


2Ie JS 




The structure of Jacobian is not of block-tridiagonal 

form as in the case of Naphthali-Sandholm method because the 

dependent variables v, . and 1, . involved in equilibrium 

Iri 3- 

relations and enthalpy balance equations are dependent on 
temperature and total liquid and vapor flow rates of all the 
stages * The corrected values of the variables are 


k +1 — k 

X = X 


+ !5 


,.k 


- - k 
A X 


(4.7) 


where k and k+1 are iteration numbers and is a damping 
parameter and X is obtained from equation (4.6). 

The following step by step procedure is employed to 
get the solution 

atep 1 : For assumed values of calculate 

F ° = F (X °) 

Step 2 1 Define inverse of the Jacobian, as matrix H which is 
updated in each iteration. As a first approximation, the 
elements of H ° are taJeen as. 
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H 

Jq is the matrix of the partial derivatives at the initial 
values of the variables. The partial derivatives are computed 
numerically as 

+ hj) - (x.) 

-■ - 

Ej 

where hj is a small increment which in the present case has 
been tal'cen to be 0,1, 


Step 3 i From the values of H ° computed in step 2 and F ^ 
computed in step 1, the correction vector ° is calculated 

as 


AX ^ = H ° F ° 


:ep 4; A value of 0 (k = 0,1,2,*.., 
Euclidean norm of 


is found such that the 


mm ***lKl. 

F(x + 0 AX ) is less than that of F (x ) i.e* 

i 

where n is the total number of equations. 

k 

Finding the value of is an iterative process in itself 

k 

and can be started with an initial value of 1,0 i.e, = 1,0 

where subscript 1 is the iteration number for finding the 
k 

optimm 0 for the Kth iteration. For the calculation of 
Ao / 0^ / .... etc, the following relations are used. 
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= [a + - 1 ] / 3 rj 

^ ''i 

H F? (X *=) 
i ^ 

This procedure is continued until the norm is redxxied. 

If after a specified number of iterations^ the norm is not 
reduced^ fresh partial derivatives are calculated at this pointy 
as indicated in step 2. 

Step 5 ; The next set of values for X is computed as follov/s; 

X = X ^ +/3 ^ A X ^ k = 0,l,2, .. .. 

The functions are evaluated with these new values of the 
variables and convergence is checked* If convergence criterion 
is not satisfied then H is updated by the lose of the following 
formulae t 

- k+1 k+1. 

5 X ^ 5 k+1 _ J k 

The convergence criterion used in this case is 

^ Pi (X) ^ £ 

i 

where €. is a jareassigned tolerance limit. 
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Step 6 s The new set of correction of variables is computed 
by following relation? 

and the whole process is iterated from Step 2 until convergence 
is achieved* 

4,1 «3 An Example 

A problem of operating distillation column processing 
mixture of eleven components has been solved using the model 
developed in previous subsection. This problem has been 
taken from Holland (1981) where solution by Thiele— Geddes method 
with © - convergence is presented, 

As seen in Table 2^ the results obtained using present 
method compare well with those reported by Holland^ thus 
establishing the validity of present method for solving multi- 
component distillation problem* 

4,2 Modelling of Crude Oil Distillation Column ; 

In the previous section a model for a simple distillation 
column was considered. In this section, a model of a crude 
oil distillation column is described* The distillation colxjmn 
with three side strippers, two pump arounds and a total 
condenser is shown in Figure 5, unlike the simple distillation 
colimn, in a crude oil distillation column there are different 
types of plates • 
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TABLE 1: STATEMENT OF THE PROBLEM 


Component 

FX. 

Specifications 

5 

2,0 

Reflux ratio ” 2.0; boiling 

2 6 

10.0 

point liqirLd feed; partial 

C3«6 

6.0 

condenser; column pressure 


12.5 

= 300 psia; no. of stages 


3.5 

= 12 and feed plate location = ! 


15.0 

reboiler temperature “ 826 °R. 

”*^5^12 

15,2 

Equilibriuim and enthalpy data 

^-^6^14 

11.3 

for all components taken from 


9.0 

Holland. 0^^*^ 

n-Cg^ie 

8.5 


400 

7.0 



"^Commonly referred to as 400^F - normal boiling fraction. 



TABLE 2: COMPARISON OF RESULTS OBTAINED IN THE PRESENT 


42 





43 


Total condenser 45® C 



FI6.5 CRUDE OIL DISTILLATION TOWER SIMULATED IN 


THE PRESENT WORK 



















(i) Plate representing condenser, 

(ii) Simple plate having no side stream and no ptimp 
around reflux. 

(iii) Plates from which liquid side stream is withdrawn 

(iv) Plates having pumparound reflux stream returned 
to them * 

(v) Plates where stripped out lighter components from 
the side-strippers are returned, 

(vi) Last plates of the colunn/side stripper where there 
is no vapor input from the plate below* 

(vii) First plates of the side strippers where feed is the 
liquid side stream draw-off from the main column. 

Besides these differences/ there is also a component (usually 
steam) in the system which is assumed to be present only in 
vapor phase and is designated as (M+l)th component* 

4*2,1 The Model: 

Each of the above types of plates will have different 
discrepency functions for material and energy balance equations. 
In this subsection, all the equations describing the plates of 
a crude distillation column are presented. 

Equilibrium Relations: 

Combining equations (2,2) and (2*4) and talcing into 
account of the fact that there is no water present in the 
liquid phase (i*e* Ij^ O# j 2/...*N) the equilibrium 

discrepency becomes. 



M 


'1 






j/i 


1 ) V . . ) - V . ^ - , 

j / i 3/ M+1 _J 


j = . 2^ * . . ,N 

If the condenser used is a total one the equilibriurn diserepency 
for first plate ( j=l) is replaced by a specification equation. 


or 


^1 “ ^1 ^■^l 

Si = V3_ ^ D 


Overall material balance equations* 
For type (1) plate 


(4.8) 


M 


j 


- h 


Vj=0 


i = 1 


For type (ii) plates 

For type (iv) and (v) plates 

Mj = L .^1 “ ° 


j = 2,5,8,li>.-12 


; j^3,6,9 


where subscript m refers to the number of sides tripper plate 
from which vapor is returned to the main column and is the 

pumparound reflux stream. 

For type (iii) plates 


^3 == ^ 3-1 


L,-V. - W.-W‘ + = O 

J j 3 3 3+1 


3 = 4/7, 10 


where Wj is the liqviid side stream withdrawn for jth plate. 
For type (vi) plates 


Mj 


L . - - L.-V. + - O 

3-1 3* 3 3 


j = 13, 15, 17, 19 


where is the steam input to the column at jth stage. 


For type (vii> plates 


= Wp L.^ - V^. + » 0 i j = 14, 16, 18 (4.9) 

where subscript p refers to the plate from which liquid 
side stream is withdrawn to be fed to the sidestrippers* 
Conponent material balance, equations; 


V. 1 , - V. . « o 

ji-#i 


7 «L 'ft • • ft ♦M 

and 3=1 


^j-l,i *“ " ^j,i ^j,i ^ ° ^ 

and j = 2,5,8,11,12 

W,+W» 

and j = 4,7,10 

W» 

l>l,i - - ''3,1 ''3+1,1 ^3+1,1 ^ = ° 

; i = 1,.....*,M 

and j = 3,6,9 
W 

Vi ■ ^J.i " '"j.i ’'j+1,1 " ° ^ ° “ 

and j =14, 16, 18 
i = 1, ,M 

and j = 13i 15, 17, 19 

(4*10) 


lj-.l,i ”* ^J,i " ^j,i - ° ^ 


Enthalpy balance equations: 

<^3-l,i ■ ’■j-i 

1=1 


ftft* V H *4" \r 

i jyi ""j+l,i j+l,i 


+ f > , h. ,) - V. H. H... = O 

3 ,i 3,1 3 ,M +1 3 ,m+l 3 + 1 ,M +1 3 + 1 ,M +1 

; j = 2,5,8,11,12 
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= 


M 

S ''j,! • ''j,i "j,i ■" ^j+i,i 

JL**** jL 

wr ’ 

+ Vi %i^ “ ^jiM+1 


’^j+l,M+l ^j+l^M+1 Vm+1 Vm+1 “ ° 

I j = 3,6^9 

where Qj is the heat removal in pumparound stream at jth plate. 


M 


'j 


^ •? 1 < 1 -s ”V. J H. . “ (1+ 

J^i/i J*”i#i i/i j#i 3 


w. + w* 


-i) 1. . h. , 


i*l 


■j 


j^i j/i 


■*■ "^j+l^i ^j+l,i^ “ Vm+ 1 ^j/M+1 ^j+l,M+l ^j+l,M+l *“ ° 


N 


/ j = 4,7, 10 


W 


Ej - V C(t^) i 4 h . - 1. . h. , - V. 

J /C^, L P/ r p#i j/i J, 1 3, 

i=l -t" 


4 H. , 
1 3#3- 


■*■ "^j+l/i ^3+1/ i^ ~ ’^jtl,M+l ^j+l,M+l " ° 


M 


Ej *= ^^3-1 i ^j-l i “ •* ^-’ -• “ '^- '• 


; j = 14, 16, 18 

. h , . - V . . H . . , 

J/i 3/3* 3/3. 3,i 


i=l 


„ tt , ^steam ^team 

“ ^j,M+l Vm+1 ^ r j = 13,15,17, 19 


^3 ‘^ta«€i^LrT\ 

where is the enthalpy of the steam being introduced 

into the 3 *th plate. 

For a total condenser (j*l) the enthalpy balance equation is 
replaced by one specification equation such as. 


Ej = - T* = O 
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where T^* is the specified condenser temperature# 

There can be other types of specifications also possible 
such as condenser heat duty# 


4#2*2 Method of Solutions 

The chosen equations (4*8, 4.9, 4.11) are solved 

simultaneously using T., L. and V. as independent variables 

3 J 3 

in the same way as for simple distillation discussed earlier 
in Section 4.1.2. The three major steps in the solution 
procedure discussed in that section are still valid except 
that in step 2, the Thcroas algorithm needs to be modified 
since the coefficient matrix structure in component material 
balance equations is no longer tridiagonal. The following 
procedure based on modification due to Hoeffling and Seader 
(1978) can be substituted for that step. Substituting the 
expression for 1 . . from equation (4.5) in equation set (4,10) 
the following relation is obtained. 


^ Vi = G . 


( 4 . 12 ) 


vhere 


''i = ''2, 

[O 


i * 




o 


"*^11,1 ~ ^12,i O ••••°] 


The structure of p matrix is shown in Table 3. 

are the elements of the tridiagonal band as usual and 

A. . and C. - . are the off-diagonal elements in this coefficient 
i/ 3 i# 3 

matrix. As shown in this table, due to the presence of 
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interconnections in the column^, there are six off-diagonal 
elements present in the matrix other than the tridiagonal 
band for this problem* Ihe equation set (4.12) is solved 
using Hoeff ling-Seader (1978) method of modified Thomas algorithm 
in which elements are forward substituted to give a matrix 
having only upper triangular elements. After fo3cward 
substitution, equation (4,12) becomes 


~ ^i (4.13) 

where /S'* is shown in Table 4, As seen from this table, because 
of the presence of off-diaganal elements, fill«!-ins occur in 
three discrete band. In this matrix main diagonal elements 
are all unity, PJs are the elements in on-diagonal positions 

ki, 

and the three of f-diagbnal bands are represented by the 
elements P. .'s. The elements of P and q. will be obtained 

1, j 1 

through the procedure described below: 



row 

1: 


row 

2: 

Bi^ (Bi ” \^i-l^ 

^i ^i ^i-l^^B^ 

row 

3: 

Same as row 2 



^3,14— S, 14-^3 

row 

4-5: 

Same as row 2 
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row 6: 


row 7-^j 


row 9 : 


rows 10-12 


rows 13: 


row 14* 


Set 


Set 


Set 


Same as rows 4—5 


Add p 


•6^16 ^ 
Same as rows 4-5 


^6^.16^®6 


Add p. , 
1/ 16 




Same as rows 7-8, 


Add p. ^ 

^9yl8 

* Same as rows 7-8 
Add p 


S, ie/®9 


” *i ^1-1/ 

Pi ^ 

as rows 10-12 


14 , i+i 


■^14.i ^1 
®14 '®14 -^4,i Pi, 14 

%4 — -V, 


'""m - Vi ®1,16 

14^ 9 9, 18 

=14 — 


) f i = 4,,.. ..12 


) 


; i = 7,, 


* • « • # 


^13 


B 


14 

^14' 

‘14 


14 


^14 

^^14 “ •^14,1 
(®14 * ■^14' 13 Pl3> 
‘^14''®14 
®14'^®14 


) 7 i = 10,... .13 


^ ^ “ 4, ...13 


14^16- 
Vi8^®14/®14 
^14 — qi4/®14 
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row 15: 


row 16 

Set 


Set 

Set 


row 17 ! 


row 18: 

Set 






18 


' >i Pi_i) 


D. 


^I6^i+1 
'16, 14 


A. 


B 


B 


16, 14 
^l6,i+l 
16 
16 


U6 


®16 

^16 

%6 

^16 
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16 




D 
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•^6,1 
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"^16 “ ^16, 13 ^13 
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^ 14 
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' ‘^16,1 ^i,l8^ 
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16 
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^16/^6 

^i-i) 

^1-1, 18/®i 
'°1 ■ \ %-l)/Bi 


^ ~ • • * <►15 
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'^a.i 
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i * lO , .,.,12 


i * 11,12 
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*18,14 <“ls'*18,:13 '’ 13 ’ 




"A-.0 • 

18 , 1 X 


i = 14 


Set 

Eis ^ 

”^18,10 ^10,16 



^18 

®18 “ •^18, i ^i,l6 

X “ 11, ••••,14 


■^18,16 

^^18 " -^18,15 ^15^ 



^I8,i+1 ^18, i ^i 

i - 16 


Set 


row 19: 


B 


B 


18 

18 


^18 

‘318 

^18 

‘^IS 

■ 

q,- ■ 


(Big - -^18,1 18^ 

^®18 " '^18^17 ^17^ 

q^Ls/^is 

‘^18 

^18 “ “^18,1 
'^IS^^IS 




10^ ...le 


10 ^ * « • • 17 


Following the above algorithm, after forward substitution, 
the matrix contains elements in upper triangle only* Applying 
backward substitution, then variables are solved as follows: 


V 


19 


^19 

- % - ^i ^i+1 
% ~ ^i ^i+1 - %18 ^18 


i = 18,17 
i = 16, 15 


^i’^ “ ^i "^i+l %18 "^18 “ ^i,16 ^16 ^ 

v^ ^ % ■" ^i '^i+1 * ^i,18 '^18 “ ^i,16 '^16 “ ’^i,14 ^14 


i = 12, 


./9 


‘^i- Vitl - %16 ^16 “ "^ 1,14 ^14 i - 8,....6 
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Vi ^ . 

''l+l - 14 ''14 

X 5, , 

Vi 

% - ''l+l 

i = 2, 1 


After having solved the component flow rates, Broyden's 
procedure can be lised as before (Section 4*1.2) to obtain 
the solution of the model equations. 

4*2,3 Flash Calculations; 

The model developed in section 4,2,1 and the method of 
solution presented in section 4.2,2 assijme that the vapor and 
liquid flow rates, temperature and compositions of feed stream 
are known. However, these properties are seldom measured 
since, before entering the fractionator, the crude oil is 
heated in a furnace and is allowed to flash on the feed plate. 

The total mass flow rate of oil, its temperature and pressure- 
at the outlet of the furnace are the only quantities which are 
known. During its course through the heating furnace the 
pipeline carrying crude oil receives heat which partially 
vaporizes the oil* This mixture flashes on entry into the feed 
plate and in order to find the temperature of feed and individual 
phase flow rates, flash calcuilations are performed according 
to following procedure. 

The temperature and pressure of the feed coming out of 
the furnace at point 1 (figure 6) is specified* Feed mixture 
in the pipeline will have both phases present at pressure of 
and temperature T^* The fraction and compositions of the 
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components in both phases in the pipeline can be evaluated 
by isothermal flash calculation as described below. When 
the temperature and pressure of the flashing zone is specified, 
the flash is termed as isothermal flash, The material balance 
and equilibrium relations of this flash calculation is given 
below: 

FXp^ = y^ (4.14) 

where E = total molar feed rate 

Xp — mole fraction of components in feed 
■''i 

Lp = liquid portion of the feed in the pipe line 
Vp = vapor portion of the feed in the pipe line 
Xj_ = mole fraction in liquid 
y^ = mole fraction in vapor. 

1 refers to the position as shown in the figure, 

= distribution coefficient 
M M 

i*=l i=l 

Equations 4*14, 4.15 and 4,16 can, be combined to give 


U - K.) 


(4.17) 


Kt + (1 - -I-) 


The equation 4.17 is solved for -|r- using Newton's method. 

If is found, x^ and y^ can be solved by using the follo^^/inc 

-p ' 1 
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equation which is obtained by combining equation 4*14 and 
4*15. 


^i = 


^ t (1 - 


{4.18) 


1 1 

substituting this value of x. in equation (4.15) y. is obtained. 
Having known the liquid and vapor fraction of feed and their 
compositions at point 1 , the enthalpies at both phases can 
be calculated at this point, fiince^ no heat is supplied from 
external agencies between point 1 and 2 the following heat 


balance equation holds good. 

F F F F 


, 2,2 2 2 

= h„ + V„ HL 

F F F F 


(4.19) 


where h - liquid phase enthalpy 
w 

Hp *= vapor phase enthalpy 

Superscripts 1 and 2 stand for locations 1 and 2 in Figure 6^ 

2 2 

Now Lp and can be found in the following way. A temperature 

T 2 is assumed and at that temperature performing isothermal 

2 2 

flash calculation as described above v_ and L is computed. 

Ihen it is checked whether equation (4,19) is satisfied. If 
not^ then a new temperature is calculated using Newton *s 
method^ and the calculations are repeated until the convergence 


is achieved 
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4,3 Thermodynamics t 

Thermodynamic properties are required to be evaluated 
for each component at different tonperatures and pressures I 

during the distillation calculations. For crude oil distillation : 

also^ it is needed to evaluate these properties for the 

pseudocomponents which represent the composition in present 

case. 

4.3.1 Correlations for Molecular Properties: 

Molecular weighty critical temperature and pressure and 
acen tricity factor of each pseudo-component are determined from 1 
the correlations given by Lee and Kesler <1975, 1976). 

Correlations ate as follows: 

MW = -12272.6 + 9486.4 SG + (4.6523 - 3.3287 SG) T^ 

+ (1-0.77084 SG - 0.02058 SG^) x (1.3437-720.79/T^) 

7 2 

10 /Tj^ + (1-0.80882 SG + 0.02226 SG ) x ( 1.8828 - 

181.98/T3^) (4.20) 

where MW = molecular weight 

SG = specific gravity, 60 ®f/60®f 

= average boiling point, “R 

= 341.7 + 811 SG + (0,4244 + 0,1174 SG) ' 

+ (0.4669 - 3.2623 SG) (4.21) 

where T ^ critical temperature, ®R 

InP^ = 8.3634 - 0.0566/SG - (0.24244 + 2.2898/SG + 0.11857/SG^ )10^'^b 
+ (1.4685 + 3.648/SG + 0.47227/SG^)10"''^T^ - (0.42019 + 
1.6977/SG )10 

where P^ = critical pressure, psi 


( 4 . 22 ) 
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w = -7.904 + 0.1352 K - 0.007465 + 8.359 

+ (1.408 - 0.01063 (4.23) 

for 0.8 

w = (In - 5.92714 + 6,09648/T^ + 1.28862 In 

- 0.169347 T^j^)/(15.2518 - 15.6875/T^^ - 13.4721 In 
+ 0.43577 T^) (4.24) 

for < 0.8 

where w = acentric factor 

K = Watson characterisation factor 

"^br * Reduced boiling point 

s 

Pj^ = Reduced vapor pressure 

4.3.2 Vapor Liquid Equilibria: 

Distribution coefficients or K-values were calculated 
assuming the system to be ideal, for which the following relation 
holds good. 

K. = pf/P (4.25) 

M- |4. 

Inhere *= distribution coefficient of component i 

s 

p. = vapor pressure of component x 

jL* 

P » total press\r:e. 

The vapor pressure data required ficr deteimiination of K- 
values are based on Maxwell— Bonnell charts in API Technical 
Data Book (1977). The following equation has been obtained 
by curve fitting and has been used in the present work. 
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following relationship: 

- Tj^ - = 2*5f (K-12) log 

where f = corrector factor. 


14.7 


(4.27) 


For normal boiling points less than 200®F (i.e. Tj^<^ 200°F) 

V'*-|>(rr ^yco^wre^ 

f=0* For all subatmosphericj^and for all substances having 

normal boiling points greater than 400®F (i.e. 400 °f 

s 

and (P /14.7) 1.0) f=l. For superatmospheric vapor pressures 

of substances having normal boiling points between 200®F and 
400°F (i.e. (P^/14.7)>1.0 and 200 °F < < 400 ®f), f is 


given approximately by the expression 

T, - 659.7 
£ = ^ 

5oo 


(4.28) 


'The ^ T is then subtracted from the true normal boiling 

1 

point, T^/ to get a new corrected normal boiling point Tj^. Nov/ 

1 s 

this new value of is used in a new calculation of P and the 

whole calculation procedure is repeated iteratively until the 

difference between two successive values of the vapor pressure 

agrees with the desired tolerance. 


4.3.3 Enthalpy: 

Enthalpy of each pseudocomponent was computed from Lee 
and Kesler Modification (Lee and Kesler: 1975, 1976; API Data 
Book: 1977) of john-Qreyson (1961) charts for evaluation of 
enthalpy. 

Liquid phase enthalpy for a particular petroleun fraction 
is given by the following relation; 
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h « A^(T-259*7) + A2(T^-259,7^) + A 3 (T^ - 259.7^) 

(4.29) 

where h » enthalpy of liqiiid petroleum fraction^ Btu/lb 
A, - 10**^ [^-1171,26 + (23.722 + 24.907 SG> + 

SG 

A 2 = 10 "^ [ (1.0 + 0,82463K) (56.086 - — ■ ■" ) 1 

A 3 = - 10 “^ 0.82463 K) (9.6757 - — 

Vapor phase enthalpy for a petroleum fraction is obtained from 
the following relations 

H 5 = h* + B- (T-0.8,T ) + B (T^-.0.64 T^) + B (T^-0.512 T^) 

■ c ^ C 3 C 

RT r — o -v- 

mr +, 5.266 w- ( S . ^y ff ) j (4.30) 

C 

where H *= enthalpy of vapor, BtiVlb 

h* « liquid enthalpy at a reduced temperature of 0,8 
calculated by eq.(4.29) 

®1 “ 10“*^ [-356 .44 + 29.72k + 84(295.02 - j 

B 2 * 10 “^ [-146.24 + (77.62 - 2. 772k)K-B4 (301,42 - J 

B 3 « 10*’^ [-56,487 - 2.95 B^} 

B 4 = - 1.0) (1.0 - (SG - 0.885) (SG-0.7){10'^)j^ 

for 10.0 < K<12,8 and 0.70 < SG < 0.885 
«= O for all other values of K and SG 
R *= universal gas constant, 1.986 Btu/(lb.mole) {®R) 

H mmt-T 

(— ~~) ** dimensionless pressure effect on enthalpy which can be 

RT^ 

found by using the following procedure. 
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< rtJ < -"IJTfr) 


( 


- H ^ (O ) I 

) mm V " o"m ' " '' V 


RT 


c,h 


‘ 3 

(4*31) 


I? (o ) 

where (g srr" ) ' *= effect of pressxire on enthalpy for the simple 

c^O 

fluid* 

,H°-H 

V Tsm .. / 




w 


®= effect of pressure on enthalpy for the 
heavy reference fluid* 

= critical temperature of the fluid for which 
enthalpy is desired, ®R 

«* acentricigiir factor of the fluid for which the 
pressvjre effect is sou^t 
= acentriciiif factor of the heavy reference fluid 
= 0*3978, 

Ihe dimensionless effects of pressure on the enthalpies of the 
simple and heavy reference fluid are to be calculated from the 


w 


(h) 


following equations 


= -1-^ 


b, + 2b3/T^ + 3b/T^ 


Ga-SCs/TC' 

r r 


T V 
r r 


d^ 


+ III + 3E 

5 T_ 
r r 


] (4. 


32) 


where 


E 


hi 


y 


P+l*“ (p+l + — ) exp (— . — ■ y ) 


r r 


Ihe superscript i refers to o for simple fluid and to h for the 

heavy reference fluid, 

T reduced ta:nper ature, T/t 
3r c 

V ~ P V/RT 
r c c 
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= ccmpressibility factor of simple fluid or 

heavy reference fluid (Z^^^). 

^ 2 ^ ^3^ ^4^ ^2^ ^3'^ ^4^ ^2^ ^ and are two sets of 
constants, one set for simple fluid and other for heavy 
reference fluid given by Lee and Kesler. 

In the present work^ pressure effect on enthalpy has 
been neglected in the first approximation. 

4*4 user^s Mannhal : 

^ computer program has been written in FORTRAN-10 and 
executed on DEC 1090 system based on the model described earlier. 
The whole program is divided and kept in six files namely MA FOR, 
FPREP FOR, XOl.FOR, PR.FOR, GEN2.F0R, FUNl.FOR. Ihe main 
program and all the subroutines used are discussed below# For 
details of the program and its usage the program listing given 
in Appendix B must be referred. 

File: MA.FOR 

This file contains only the main executive program. 

All the necessary inputs have to be supplied in this file either 
through input data file (F0R20,DAT) or data statanents. Initial 
guesses of the variables are also read through sn input data 
file (F0R40.DAT). On successful convergence, it gives the 
values of the variables at solution which can be either stored 
in an output file (OUTiDAT) or can be displayed on the 
terminal screen. The properties of the pseudo-qomponents 
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are also printed or mu. . 

>^®Pl^yed. Ihe Input variaMos for this 
PJaogram are as follows: 


IPRINT ! 


NO 


II^AX; 

M: 

N: 

P: 

TBPs 

API: 

WPR: 

FEED: 

TFEED; 

PPEED: 

NFEED: 

DIS: 

REFLUX: 

WL: 


Printing control parameter 

= 1.- When printing of input data is needed, 
number other than 1 win skip printing these dal 
Output control parameter 

= 05,- When results are needed to be displayed on 
the screen of a terminal 

= 22r When results are required to be stored in 
file OUT.DATs 

Maximum number of iters +-i- ram <=. 

Iterations permissible. 

No, of components 

No, of plates 

Operating pressure, psia 

Average boiling point of each dsgu,^o 

'-d'-n pseudo component, ®c 

Gravity of each pseudo component 
height fraction of each pseudo component 
Total mass flow rate of feed, kg/ti 
Temperature of the feed, “c 
Pressure of the feed, psia 
Feed plate number 
Distillate flow rate, kg/h 
Reflux ratio 

Liquid side stream draw off, kgmol/h 
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WLI Sc WL2: Flow rates of pumparound reflux streams, kgmol/hi 
Ql & Q2s Heat ranoval from pumparound reflux streams, Kcal/h 
STM: Rate of steam input to the column/sidestrippers, Kg/h 

File: FPREP.FOR 

This file contains subroutines which are needed for feed 
characterization* 

Subroutine FPREP: 

This subroutine computes the molecular properties of 
pseudocomponents and evaluates the vapor and liquid fractions 
of feed and their composition through flash calculation* 

Input data: 

FEED, T, P, M, WFR, API, TBP 
Where T: temperature of feed, °C 
P: pressure of feed, psia 

Output data: PC, TC, W, AMW, SG, AMWFD, FL, FV, FLC, FVC, ENTHL, 
ENTHV 

where PC: Critical pressure of each pseudocomponent, psia 
TC: critical temperature of each pseudocomponent, ®R 

W: acentriciHy factor 
WK; Watson characterization factor 
AMW: Molecular weight of each pseudo component 
SG: Specific gravity 60®P/60‘’F 
AMWFD:Average molecular weight of feed* 

FL: Liquid portion of feed, kgmole/h 
FV; Vapor portion of feed, kgmole/h 
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FLC: Composition of liquid phase of feed> kgmole/h 
FVCj Composition of vapor phase of f eed^ kgmole/h 
ENTHL: Enthalpy of liquid phase of feed, kcal/h 
ENTHV: Enthalpy of vapor phase of feed, kcal/h 
Subroutine FLASH: 

This subroutine is used for isothermal flash calculation 
Input Data: FM, XF, T, P, M 

where FM: molar flow rate of feed, kgmol/hr 
XF: mole fraction of feed 

T: temperature at which flash calculation is sought, °C 
P: pressure of flash zone, psia 
Subroutine EQMl 

This subroutine computes ideal K-values using Maxwell- 
Bonnell vapor pressure correlations* 

Input data: T, P, WK, TBP, M 
Output data: QK 

where QK : K-value of each component 
Subroutine FLl 

This subroutine performs the adiabatic flash calculation 
Input data: FM, XF, P, M 

Output data: ENTHL, ENTHV 
File : XOl^FOR 

File X01*F0R contains all the subroutines necessary tc. 
solve component material balance equations* 
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Subroutine XOl 

This subroutine generates the coefficient matrix of 
component material balance, equations and then solves calling 
modified Thomas algorithm subroutine (THOi^iAS) 

Input data : T, N, M, Py REFLUX, RL, V, WLl, WL2, WL, PDC, PVC 

where T; temperature of each stage, °C 
N : total number of stages 
RL: liquid flow rate at each stage, kgmol/h 
V; vapor flow rate at each stage, kgmol/h 
Output data: RLC, RVC 

where RLC : component flow rate of liquid at each 
Stage, kgmoleAi 

RVC: component flow rate of vapor at each stage, 
kgmoleAi 

Subroutine THOMAS 

Input data: AT, BT, CT, F, OFl, 0F2, 0F3 

v;herG 7».T, BT, CT: elements of tridiagonal band of 
the coefficient matrix 
OFl, 0 f 2, 0F3: Off diagonal elements in the 
coefficient matrix 

F : right hand side of the component material 
balance equation 

File: GEN2.F0R 

This file contains the following subroutines. 
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Sxabroutine GENG: 

Subroutine GENG computes the discrepancy function of 
jth plate* 

Input data: J/ NJ, XK^M, REFLUX/ FL, FV, ENTLF^ EMTVF, WL, 
WV^ STEAM, WLl, WL2, TCCN'4 RLC, RVC 
where J : jth plate 

N : total number of equations per stage i.e. 3 
NJ : number of stages 
XX : vector of independent variables 
Output data; G : 

where G i the vector of discrepency functions at jth 
stage . 

Subroutine BLOCK 

Subroutine BLOCK returns the values of independent 
variables at solution on successful exit. 

Input data; XX, N, NJ, M, P 

whore Kj: number of stages 

N ; number of equations per stage 
Output data: XX 

where XX ; the vector of independent variables 
Subroutine JAC 

This subroutine evaluates the partial derivatives of tl 
Jacobian matrix using numerical differentiation* 

Input data; XX, Gl, N, M, P 



71 


where XX s the values of the independent variables at which 
partial derivatives are sought, 

Gl : values of the discrepencies at that point 
N : number of equations per stage 
Output data: PD 

where PD: partial derivatives 
Subroutine Yecmut 

This subroutine performs the multiplication of a matrix 
with a vector. 

Input data; X< Y, N 

where X : matrix 
Y ; vector 

N ; order of X matrix 
Output data : Z 

where Z : resultant vector 
File: FUNl.FOR 

The file contains only one subroutine, 

Si±)routino FlU 

This subroutine computes the discrepencies of all the 

stages. 

Input data: XX/ N/ M, P# I'TL/ WLl/ WL2 

where N: total number of stages 

XX: vector of independent variables 
Output data : GT 

where GT: the matrix whose elements contains the 
discrepancies of all stages. 
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File: PR .FOR 

File PR. FOR contains the all the subroutines required 
for thermodynamic property estimation i.e, vapor*<-liquid 
equilibria and enthalpy. 

Subroutine EQILl: 

Ihe subroutine computes K- values (Ideal) using Maxwell- 
Bonnell correlations for vapor-pressure. 

Input data: J, T, P, WK, TBP 

vAiere J : jth stage 

T : temperature of jth stage# 

P ; pressirre of jth stage# psia 
Output data: RK 

where RK : K- values of each component at jth stage* 
Subroutine ENTHLP 

This subroutine computes the enthalpies of pseudo- 
components for both phases. 

Input data : J# T# P, M# WK# SG 
where J : jth stage 

T : temperature of jth stage, °C 
P : pressure of jth stage# psia 
Output data: HL# HV# CPL# CPV 

where HL : liquid phase enthalpy of each component at 
jth stage# kcal/kg mol 

HV ; vapor phase enthalpy of each component at j 
stage# kcal/kg mol. 



CPL : specific heat of each component in liquid 
phase at jth stage, kcal/kg mol°c 
CPV : specific heat of each component in vapor phase 
at jth stage, kcal/kg mol°C 
Subroutine SHEAT 

This subroutine computes enthalpy of steam 
Input data: J, T 

where J : jth stage 

T : temperature of jth stage, ®C 
Output data: HSTM 

where HSTM ; enthalpy of steam at jth stage, kcal/kg mol. 



CHAPTER 5 


RESULTS AND DISCUSSIONS 

Using the model developed in the previous chapter/ an 
atmospheric crude oil distillation column has beem simulated 
and the results have been compared with those obtained from a 
commercial simulation package* Also probable reasons are 
discussed for the discrepancies observed* 

5*1 Statement of the Problem ; 

The column configuration has been shown in Figure 5. 

The main column has 13 equilibrium stages (plates) and the three 
side strippers have 2 plates each. At the bottom^ of the main 
column and each of the side strippers^ open steam is introduced 
at 3.4 atm pressure and 350‘’C. Steam flow rates are given 
in the figure. Feed at a rate of 100 tonnes/hr at temperature 
340°C and pressure 14 atm is introduced at plate 12. The light 
ends analysis of the feed is given in Table 5. The TBP and 

specific gravity data of the crude is given in Table 6. The 

column is operated at 4 atm press uire. The main coltann has a 
total condenser/ the temperature of which is 45 “C. Distillate 
flow rate is 160 kgmole/hr* There are two pumparotind redEluxes 
employed in the main column* From the 7th tray of the tower 

a liquid fraction of 616 kgmole/hr is withdrawn and returned 

to 6th plate after extracting 2.75 x 10^ Kcal/hr of heat from it. 
In the second pumparound reflux, a liquid portion of 275 kgmole/hr 
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is withdrawn from 10th plate and returned to 9th plate after 
cooling it, Ihe heat removal from this reflux is 2,0 x lO^Kcal/hr. 
A liquid sidestream of 98,43 kgmole/hr is withdrawn from 4th 
plate and fed to the top plate of first sidestripper and the 
vapor leaving the top plate of this side stripper is led to 3rd 
plate in the main colunn, 178,84 kgmoleAir of liquid withdrawn 
from the 7th tray of the main colunn is the feed to the 2nd 
sidestripper and vapor coming out of this stripper is introduced 
to plate 6 of the main column* The feed to the last sidestripper 
is a portion of liquid withdrawn from plate 10, The flow rate 
of this withdrawn stream is 129,15 kgmole/hr, Ihe stripped vapor 
from 3rd stripper is sent back to the 9th tray of the main 
column, 

5*2 Results and Discussion : 

Molecular properties of the pseudo componeribs as 
calculated by the correlations of section 4,3 are given in Table 7* 
In this table the components have been characterized by their 
average boiling points. The liquid fraction and vapor fraction 
of the feed and the composition in both phases as computed from 
flash calculations are tabulated in Table 8, In the same table^ 
the results from the standard solution are also included for 
comparison. In view of the mismatch, the flash temperature vj’as 
increased by 9 ®C which brought the total vapor and liquid 
flow rates close to the standard solution. The vapor flow rate 
was 468,96 (470*15 from standard solution) and the liquid 



flow rate was 94.51 (96.04 from standard solution). Two sots 
of calculations were made, (a.) with feed to the column provided 
from the flash calculation with increased tenperaturQ (b) with 
feed same as that obtained in standard solution. 

Table 9 shows a comparison between calculated plate 
temperatures and total vapor flow rates obtained for the above 
two cases and the standard solution, as seen from this table 
deviations from standard solution in case (a) are larger in 
comparison to case (b). This is to be expected since case (a) 
represents the added discrepancies of both flash calculations 
and distillation calculations, whereas case (b) shows only 
distillation column calciolation errors. The possible reasons 
for these deviations sea:ns to be thermodynamic in nature. 

The vapor enthalpies calculated in the present program are 
higher than actual since the pressure correction was not included 
This leads to a lower vapor fraction in the flash zone as shown 
in Table 8. By increasing the temperature of the flash, the 
vapor fraction is increased to match the standard solution 
but then, the heat content of vapor will also go up thereby 
requiring higher reflux. This amounts to higher vapor flow 
rates through the column as seen in case (a) af table 9. The 
match between case (b) and the standard solution is more 
reasonable, all deviations being less than 10 percent. Here 
again, the enthalpies may be responsible for the observed error , 
The vapor -liquid equilibrium constants or K-values have been 
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calculated by application of the idealized JRaoult’s lav/ 
where their dependence on composition is not included. 

“Ibis may have also contributed to the loss of accuracy* 

% 

Tables 10 and 11 show a comparison of product compositions 
obtained in present study for cases (a) and (b) with those 
of standard solution, as expected, the match between case (a) 
and standard solution is poorer as compared to that between 
case (b) and the standard* 


I 



78 


TABLE 5 s LIOIT ENDS ANAL-^SIS 


Component 

Propane 

I-Butane 

N-*Butane 


Wei^t percent 
0.1 
0.1 

0.3 


I-Pentane 


0.3 


N-Pentane 


0.5 
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TABLE 6 : TBP AND GRAVITY (API ) DATA 


No, 


Cut ®C Cumulative wt,% 


API 


1 

Upto ; 

25 

1,9 

2 

25 

- 

50 

3.3 

3 

50 

- 

75 

7.2 

4 

75 

- 

100 

11.9 

5 

lOO 

mm 

125 

17,4 

6 

125 

mm 

150 

22.6 

7 

150 

- 

175 

29.1 

8 

175 


200 

33.1 

9 

200 

- 

225 

36.8 

10 

225 

- 

250 

41.4 

11 

250 

- 

275 

47.4 

12 

275 

- 

300 

51 .8 

13 

300 

- 

325 

56.8 

14 

325 

- 

350 

61,6 

15 

350 


370 

65,4 


95,00 

85,30 

63.60 
55.98 
52,08 
48.02 
47.28 

45.60 
41.76 
39.12 
34.70 

33.22 

36.22 
33.58 


32.42 


molecuiar properties of the pseudo-compc^ehts 
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CHAPTER 6 


CCtJCLUSia-^S AND RECOMMENDATIONS 

The mathematical model developed in the present study 
for multicomponent fractionation has been found to be 
computationally stable and effective* The model gave results 
comparable to those reported by Holland for a 11 component 
problem* Ihis method has been extended to solve crude oil 
distillation problem using pseudo— component concept and the 
results have been compared with those obtained with a standard 
package* The probable reasons for the discrepencies have been 
discussed. Pressure gradient throughout the colunn which 
was neglected for the sake of simplicity may actually influence 
the results. It is hoped that making necessary corrections 
to account for effect of pressure on enthalpy and using an 
appropriate nonideal equation of state for VLE calculation 
will improve the results, 

Ihis work can be extended on following lines to develops 
a complete simulation package for the crude oil distillation, 

i 

In the work initial guesses for all independent variables were 
supplied to solve the problem. But to make the program more 
effective^ an in-built procedure to generate the initial 
guesses should be developed. This can be done by making use of 
empirical short cut methods of crude oil distillation calculatit. 
such as one discussed by Watkins (1979), Using this method 
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temperature^ vapor and liquid flows at all key trays of the 
tower can be computed and then by linear interpolation they 
can give rise to the initial guesses for the independent 
variables at all stages. In the present work, the liquid 
streams withdrawn from the main column (kg moles/h) were 
taken as specification. But in actual practice the products from 
the side strippers (kg/li) will be the specifications. This 
can be obtained by modifying the subroutines for calculation 
of discrepancy functions# Another irrportant specifications in 
crude oil distillation is amount of overflash. This also should 
be properly accounted for • Since efficiency, in case of cmde 
oil distillation, is very difficult to define, an analogue 
colunn has to be found by trial which contains equilibrium 
stages but performs like the real column# Since Naphthali- 
Sandholm method can handle non-ideal thermodynamics, can 
incorporate efficiencies and is known to have better convergence 
characteristics, it may be worthwhile to approach the 
modelling of crude oil distillation using this method on a 
machine which has sufficient memory. 
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APPENDIX A 


NAPHTHALI-SANTHOLM FORMULATIOJ OF CRUDE OIL 
DISTILLATION PROBLEM 


For Naphthali-Sandholm method the functions to be 
solved are component material balance equations for each plate^ 
the equilibrium relationships and enthalpy balance equation 
for each plate* The material balance and enthalpy balance 
equations are same as Eqs* 4*12 and 4*11 • Only equilibrium 
relations will be described in this section, 

V 




E..K..1.. V.,.+v 

■ Jy .3 - - (1 - E, p 

j •’ j +1 m 


(A,l) 

for j = 3, 6^ 9 

where m is number of the plate other than (j+l)th plate from 
which vapor is introduced to the jth plate of the column. 


V , . 

V . 


K 




j/i h. 


j = 13, 15, 17, 19 

V . 


j/i 


4^^ - E . . K . . -1^ - ( 1- E . . ) -4— 

Vj 3,1 3,i Lj 3,1 ^3+1 


(A.2) 

(A.3) 


for all other plates. 

Since the last plates of the column (13, 15, 17, 19) do not have any 
plate below them, it is not possible to define Miorphree— plate 
efficiency for them. That is why the efficiency in this case 
has been assumed to be unity. In this method dependent variables, 
total liquid flow rate and vapor flow rates, are calculated by 
summing up the component flow rates. 
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Solution techniques 

The functions and the variables are grouped stagewise# 
They are solved simioltaneously by n-R technique where 

J ^ = - F (x^) 

and *= X^ +^‘3 Ax^ 

The Jacobian matrix has got the same structure as that shown 
in Table 3* Only difference is that, the elements of matrix 
Shown in Table 3 are algebraic numbers, where as in this ca.se, 
they will matrix blocks of order (2 m + 1). The Jacobian matri:>c 
can be solved using modified block Thomas algorithm* The 
algorithm developed in Section (4*2) to solve the Component 
material balance equations is used* In the developed algorithm 
as shown, the , operations like addition, subtraction, division 
are algebraic as the elements are algebraic numbers, but in 
this Case there Wiii ha matrix operations as the elements of 
the Jacobian are matrices* 

To facilitate the convergence, one damping factor is 
used Which is calculated using Golden Section Search technique • 
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PROGRAM LISTING 



